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Abstract 

Encouraged by the AdS/CFT correspondence, we study emergent local geometry in 
large N multi-matrix models from the perspective of a strong coupling expansion. By 
considering various solvable interacting models we show how the emergence or non- 
emergence of local geometry at strong coupling is captured by observables that ef- 
fectively measure the mass of off-diagonal excitations about a semiclassical eigenvalue 
background. We find emergent geometry at strong coupling in models where a mass 
term regulates an infrared divergence. We also show that our notion of emergent geom- 
etry can be usefully applied to fuzzy spheres. Although most of our results are analytic, 
we have found numerical input valuable in guiding and checking our results. 



1 Introduction 



The emergence of geometry has long been recognised as a key issue in quantum gravity. The 
more recent discovery of the AdS/CFT correspondence [1] has indicated a complementary 
viewpoint: Certain large N field theories are most simply described starting from a higher 
dimensional dynamical geometry. The geometry is not apparent in the weak coupling 
(Lagrangian) description of the quantum field theory, and in this sense is emergent. The 
emergence of spacetime can be thought of as a precise realisation of 't Hooft's insight [2] 
that large ./V gauge theories are string theories, together with the fact that string theory 
describes a theory of quantum geometry. 

The best understood cases of emergent geometry from field theory begin with D branes 
in a pre-existing geometry 1 . The geometry 'dual' to the strongly coupled field theory on the 
D brane worldvolume is then obtained by computing the gravitational backreaction of the 
D branes [1, 4]. These theories are often supersymmetric and, in the best understood cases, 
conformal. One would want a much more comprehensive understanding of how geometry 
emerges, and the relevant situations where it can be applied. In essence, one needs to find 
a way to solve a key aspect of the field theory dynamics and find geometry. 

There is ample reason to believe that the large N expansion might be relevant for 
QCD, the theory of the strong interactions, and one would like to have a first principles 
approach to calculating the corresponding dual geometry (or dual string theory). If this 
is understood, exploiting the geometrical information may allow us to describe the strong 
coupling dynamics of the theory in a more economical way. The strong interactions are not 
supersymmetric and, in the regime of strong coupling, they are not conformal either. 

In this paper we would like to study emergent geometry from first principles in a sim- 
plified setting. We will study large ./V systems in zero dimensions, with and without su- 
persymmetry. We will solve the multi-matrix models exactly in certain limits, and look for 
emergent geometry. 

The simplest emergence of geometry from large N matrices occurs in a Gaussian matrix 
model for a single matrix. In the large N limit, the integral is dominated by a saddle point 
in which the eigenvalues of the matrix are distributed in a semicircle [5]. This semicircle 
can be thought of as a continuum geometry that emerges at large N. This very simple 
example already illustrates an important theme for us. The emergent geometry is possible 
because the Gaussian mass term balances the repulsive inter-eigenvalue force. The Gaussian 

1 Early pre-stringy observations of emergent geometry include for instance [3]. 
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model may be generalised by introducing interaction terms for the matrix, and an elegant 
mathematical theory allows us to find the eigenvalue distribution in that case [6]. 

A tractable step beyond single matrix models are normal matrix models, in which a 
Hermitian matrix and its Hermitian conjugate commute with one another. These models 
show an emergent geometry at large N which describes a two dimensional droplet [7] . Again 
there is an elegant mathematical framework to describe these models, see for instance [8, 9]. 

For general multi-matrix models, even with just two matrices, the question of whether 
there is or not an emergent geometry is difficult and a developed framework is lacking. 
Unlike the case of single matrix or normal matrix models, these systems can rarely [10] be 
solved exactly even at large N. Some recent numerical work can be found in [11, 12, 13, 14, 
15, 16]. However, many of the cases of most interest are numerically problematic because 
of the sign problem in certain supersymmetric theories. 

The AdS/CFT correspondence [1] suggests that a natural organising principle is a strong 
coupling limit. For the correspondence to work, the degrees of freedom in the field theory 
need to 'geometrise' at strong coupling in order to reproduce the dual higher dimensional 
(dynamical) spacetime. There are various discussions of emergent geometry in AdS/CFT 
in the literature, including [17, 18, 19]. However, only recently has an attempt been made 
to systematise the emergence of geometry as a consequence of the strong coupling limit [20] . 

With multiple matrices, it seems that a key aspect of an emergent classical geometry is 
that the matrices commute with each other in the large N limit. In this way the typical 
iV 2 degrees of freedom of matrices get effectively reduced to order N degrees of freedom 
at low energies. The collective description of these low energy degrees of freedom can 
often be given in terms of a joint eigenvalue distribution for several matrices. It is the 
geometrical description of this eigenvalue distribution that produces the emergent geometry. 
One objective of this paper is to make the notion of commuting matrices more precise. In 
[20] it was proposed that the zero modes of the six scalar fields of M = 4 super Yang-Mills 
theory on a spatial S 3 commuted at strong coupling. From this proposal one can show that 
the joint eigenvalue distribution of these six matrices forms an S 5 that should be identified 
with the geometric S 5 that arises in the AdS§ x S b of the dual IIB string theory. 

The proposal of [20] was subsequently generalised to describe orbifolds of the N = 4 
theory [21] and also to M = 1 theories [22, 23]. Various successful checks of the proposal were 
performed in [24, 23, 25]. However, given that the theory involves infinitely many coupled 
matrices with comparable masses, it is difficult to prove the validity of the truncation to just 
six fields. Furthermore, the theories considered have all been superconformal and somewhat 
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similar to the maximally symmetric Yang-Mills theory in four dimensions. 

In this paper we will look at models in which we can show directly and unambiguously 
whether the matrices commute or not at strong coupling. The paper is organised as follows. 
We will solve a succession of multi-matrix models in the strong coupling limit. At each step, 
we shall check our analytic results with numerics. Although many of our final results are 
analytic, the interplay with numerics has been a very important guide towards these results. 
In a later section we consider some models that we cannot solve analytically. Finally, we 
also show how our notion of emergent geometry can be applied to fuzzy spheres. The 
concluding discussion includes a summary of our results. The main point we emphasise is 
that higher dimensional emergent geometry can arise naturally in models where a mass term 
regularises an infrared divergence of a massles model. In these cases, at strong coupling, 
the eigenvalues are sufficiently spread out that off-diagonal modes do not contribute to low 
energy dynamics. 



2 A two matrix model at strong coupling 
2.1 Solving the model: Parabolic distribution 

We begin by considering a two matrix model that can and has [26] been solved exactly at 
all couplings. We could read off most of the quantities we need from [26]. Instead we will 
solve the model in a more low-tech way, that is geared towards the more general issues we 
would like to understand at strong 't Hooft coupling. We will see that this model can be 
re-interpreted as an emergent two dimensional geometry at strong coupling. 
Consider the two Hermitian matrix model 

Z = J vXVYe- trX2 ~ trY2+ s 2 ^,Y] 2 . (!) 

Because this integral is quadratic in both X and Y, we can diagonalise X and then perform 
the Y integral exactly. The partition function becomes 

Z = fdx!... dx N e~ Zi* 2 TT — J ? ~ x i (2) 

= J dx! . . . dx N e~ ^ x2+1 2 ^ iog[i+ 9 2 (^-^) 2 ] _ (3) 

The first logarithm appearing here is the standard Vandermonde determinant arising upon 
diagonalising X. We can look for a large N saddle point to this integral. The saddle point 
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equations of motion are 

These equations have been solved at large N by Kazakov-Kostov-Nekrasov [26], somewhat 
implicitly. Here we shall be interested in the strong 't Hooft coupling limit A = g 2 N » 1, 
and shall take a less sophisticated approach. 

In the large N limit, the saddle point equation becomes 

p(y)dy / r\ 

( X -y)[l + g2( X - y y}- { > 

We took the continuum limit ~^ f p{x)dx and introduced the eigenvalue density p(x). We 
are working with the normalisation f p(x)dx = N . As usual, a principal value is understood 
in the integral (5). 

We can solve this equation at strong 't Hooft coupling by first noting the general result 

p(y)dy 7rp'{x 



■L 



(x-y)[l + J 2 (x-y) 2 } J 
One can check that indeed 

(1 - y 2 )dy 2ttx 
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+ • • • as J-too. (6) 



H as J^oo. (7) 



(x-y)[l + J 2 (x-y) 2 ] J 
The key point is that this result is linear in x. From this integral, it follows that in the 
strong 't Hooft coupling limit, A = g 2 N » 1, the solution to the integral equation (5) is 
the parabolic distribution 

, , 3iV(L 2 - x 2 ) 
P{X) = 4L3 ' (8) 

with 

2 l/3 A l/6- ^> 

The correctness of this solution can be verified by plugging it into the equation (5), perform- 
ing the integral and then taking the large A limit. We have also simulated the eigenvalue 
partition function (2) numerically using a Hybrid Monte-Carlo algorithm. The resulting 
distribution is shown for N = 500 and A = 600 in figure 1. We stored 10000 configurations 
and determined the distribution p{x) from the distribution of x 10000 points. The plot 
shows that the eigenvalue distribution is indeed close to our theoretical result (8). An anal- 
ogous plot at lower A (say A = 20) shows clear deviations away from being a parabola. Note 
that the appearing in these computations is just a discretisation of the integral equation 
(5). Therefore, for these purposes, we can take A to be larger than if we wish, without 
upsetting the 't Hooft limit. 
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Figure 1: Numerically simulated eigenvalue distribution solution to (4) with N = 500 and 
A = 600, together with the theoretical result (8), which has L = 12.91 in this case. 

We should note however, that the solution (8) is not correct very near to the endpoints 
of the distribution, x ~ ±L. More specifically, it is not correct for (L 2 — x 2 ) < L 2 /A 1 / 3 . 
This occurs because there are higher order corrections to (7) of the form 1/[(1 — x 2 )J] n , 
with n > 0, which become large near the endpoints. However, this represents a vanishingly 
small fraction of the support of the eigenvalue density and so should not be important for 
many observables. Specifically 

# incorrect 1 



0. 



(10) 



# total A 2 / 3 

Furthermore, by going to very large A we have checked that the result for the width L in 
(9) agrees excellently with numerics, including the prefactor (37T/2) 1 / 3 ~ 1.68. 
For a typical pair of eigenvalues, the solution (8) implies that 
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: (x-y) 2 ~AL 2 /N~A 2/3 >l. 



(11) 



This observation suggests that the matrices X and Y are effectively commuting at strong 
coupling, because most of the off diagonal modes are parametrically more massive than the 
diagonal modes. Recall that the off diagonal mode connecting the ith and jth eigenvalue 
has mass 1 + g 2 (xi — Xj) 2 , as we used for instance in evaluating the determinant (2). We 
shall make this notion of commutativity more precise in the following subsection. 
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It is now easy to evaluate, for instance, 

tiX 2 _ f L 3x 2 (L 2 -x 2 )) ^ _ L 2 _ (3ir) 2 / 3 N 
N ~ J_ L 4L3 dx ~ 5 " 5 • 2V3 AV3 + ' ' " ' {U > 

We can note that this agrees exactly with the strong coupling result by Kazakov-Kostov- 
Nekrasov, see [26] equation (6.29). 

2.2 Commutators and criteria for locality 

Let us now consider observables that depend on both matrices X and Y in the two matrix 
model. Our objective here is to show more quantitatively that the two matrices commute 
at strong coupling. Our comment in the previous subsection about the large mass of off- 
diagonal modes is insufficient in itself. In the framework we are using so far, in which X 
is diagonalised and Y treated exactly, the off-diagonal modes do in fact make significant 
contributions to generic observables involving Y. We shall see this shortly. When X and 
Y are treated on an unequal footing, it is best to discuss basis-independent quantities. In 
particular, we are interested in the following combinations 

tr{XYXY) , ti{X 2 Y 2 ) , tr[X, Y} 2 = 2 (tr(XYXY) - ti(X 2 Y 2 )) . (13) 

Obviously the commutator square should have information on how close to commuting 
are two sets of matrices. The other two combinations have different large N behavior in 
free matrix models. The term tx{XY XY) would vanish at the planar level, and tx{X 2 Y 2 ) 
would not. We would consider the case of completely uncorrelated matrices to mean that 
the matrices are non-commuting. The ratio 

tr(XYXY) 

tx(X 2 Y 2 ) ' { ' 

would then serve as an order parameter that tells us something about the correlation of the 
matrices. To leading order in l/N it vanishes in the Gaussian model. 

It is easy to show that \r\ < 1 for Hermitian matrices. It results from a simple manipu- 
lation of the following two inequalities: 

tv([X,Y} 2 ) <0, (15) 
tr{{X,Y} 2 ) > 0. (16) 

These follow from the fact that [X, Y] and {X, Y} are anti-Hermitian and Hermitian respec- 
tively. For matrices that commute, we find that r = 1, while for matrices that anticommute 



6 



we would find that r = — 1. In the case r = — 1, the square of the matrices would commute 
with each other, so there is also a lot of order in the eigenspaces of the matrices. 

When we consider our two matrix model, the easiest of these to compute is the commu- 
tator squared. Specifically 

Here log Z should be evaluated on the large N saddle 

\ogZ = - j p{x)x 2 dx + j dxdyp{x)p{y) Qlog(z - y) 2 - ^log[l +g 2 (x - yf 
Using our parabolic solution (8) we obtain to leading order at large A 



(18) 



TV 3 

tr[X,Y} 2 = __ + .... (19) 

This is the correct answer for the commutator to leading order in large A. It is clear from 
the computation that higher order corrections to the eigenvalue distribution give subleading 
contributions. 

To get the other traces we can introduce a source Jjj into the action for the Yy compo- 
nent of Y. Let 

Z[J] = J PIWe 4rl2 -" y2+ 9 2tr ^l 2+lrjy . (20) 
Diagonalising the X matrix, this becomes 

Z\J\ = el^i^Jii^+a 2 ^-*^- 1 Z[Q] . (21) 

Thus 

t W yv yy) _ \^ x i x J b 2 Z[3\ f dxdyp{x)p(y)xy 

tr(XYXY) - ^ -J ^ _ , (22) 

and similarly 



dxdyp{x)p(y)x l 



+ g 2 (x- y ) 2 Y (23) 

It is easy to evaluate these integrals using our large N distribution (8) to obtain at large A 

lvvvvs 3^ /3^\ 1/3 iV 3 3iV 3 aN 3 
tT(XYXY) = — (— J _-_ + _ + ..., (24) 

37r {3it\ 1/3 N 3 N 3 , aN 3 



In these expressions we have included an unknown contribution aN 3 /X that comes from the 
leading correction to the parabolic eigenvalue distribution (8) in the large A limit. We know 
that the contribution has to be equal in the two expressions because taking their difference 
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recovers our previous result for the commutator (19). Recall that our expression for the 
commutator did not depend on corrections to the eigenvalue distribution. The ratios we 
are about to consider do not depend on a. For completeness, in Appendix A we show that 

We checked the leading order results for the commutators numerically, simulating the full 
partition function (1). 

We can now ask how these observables capture the commutativity of the matrices at 
strong coupling. The most naive object to look at would be 

NtT[X,Y]* _ 25/2^ 1 

trX 2 tr y2 " 2 [sttJ A 1 / 3 + [ ' 

This says that the commutator is vanishing relative to observables that only depend on 
single matrices. The ratio cancels out the overall scaling of the X and Y matrices. Thus, 
the condition that (27) go to zero seems to be a natural notion of whether the matrices X 
and Y commute. As we discussed above, another natural ratio to consider is 

tr(XYXY) _ 35 / 2 y/ 3 1 

tr(A2y2) - 1 67T { 37T ) X W + " W 

Which also shows that the matrices commute. The parameter A 1 / 3 controls the size of the 
non-commutativity of the matrices. If the matrices are sufficiently close to commuting, we 
make small errors by assuming that they are mutually diagonal, and that there is a joint 
eigenvalue distribution. 

Therefore, we consider the behaviour of the ratios (27) and (28) as criteria for the 
emergence of a local geometry. Our computations above show that they are very closely 
related to the mass of the off diagonal modes, m^ d ~ A 2 / 3 , as we should expect. 

One should notice that the criterion for commutativity, r can be refined further. For 
example, we can ask more local questions in the spectrum of X if we consider the ratios 

_ trU(X)Yf(X)Y) 
f ~ tr(/(X) 2 Y 2 ) ' { > 

for / some real function of x which is peaked in some region. We can do the same with 
Y. Depending on how rj varies with the width 5 of /, we can talk of a local degree of 
noncommutativity on the scale of S. It should be noted that these can be easily evaluated 
numerically if / is a rational function, without the need to diagonalize X. This seems to 
serve as a reasonable definition of the local sharpness of a fuzzy geometry. 
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2.3 Emergence of local geometry: hemisphere distribution 

Given that most of the off diagonal modes of Y are parametrically heavy when A 3> 1, and 
given the observations of the previous subsection, we might expect to be able to recast this 
model as a commuting matrix model for the two matrices X and Y . A two dimensional 
commuting matrix model is a matrix model in which the matrices are further constrained 
to commute. At large N the model is described by a joint eigenvalue density p(x, y) for the 
eigenvalues of the commuting matrices. 

It is easy to see that the hemisphere distribution 



p{x,y) 



3Nt/iJ~ 



y 



o 



2vrL 3 
otherwise , 



for x 2 + y 2 < L 2 



(30) 



recovers our one dimensional parabolic distribution (8) upon integrating out one direction 

rVL 2 ~x 2 

p{x) = / p(x,y)dy. 



(31) 



This immediately implies that all observables depending on only one of the matrices can 
be computed using this two dimensional eigenvalue distribution, which we might call the 
hemisphere distribution. Note that (30) is the unique radially symmetric distribution with 
the property (31). Radial symmetry is appropriate as the original two matrix model was 
SO (2) invariant. The emergent hemisphere distribution is shown in figure 2. 




Figure 2: The emergent two dimensional hemisphere distribution. 
A nontrivial test of this emergent two dimensional eigenvalue distribution is to reproduce 



9 



observables that depend on both X and Y. In particular we find 

/3-7T /3"7r\ 1 ^ /V 3 
p(x, y)x 2 y 2 dxdy = — J ^ . (32) 

Precisely reproducing the exact result (24) to leading order at large A. Here we note that 
off-diagonal elements are manifestly not necessary in order to correctly compute observables 
involving both X and Y in this commuting framework (except for commutators, of course, 
which vanish to leading order). Unlike in the previous discussion in which we diagonalised 
only X, in this simultaneously diagonalised basis there is a well-defined separation into light 
eigenvalues and heavy off-diagonal modes. 

It is natural to ask for an action that has the eigenvalue distribution (30) as its ground 
state. The naive thing to do is to obtain a one loop effective action for the simultaneous 
eigenvalues of X and Y by integrating out the off diagonal modes. One might hope this 
will work even at strong coupling, because, as we saw, the off diagonal modes are becoming 
parametrically heavy. An estimate of when perturbation theory is valid can be obtained 
as follows (essentially this argument appears in, for instance, [15]). Given an eigenvalue 
distribution of extension L, the action for off diagonal modes SX is schematically 

SWdiag ~ (1 + g 2 L 2 )5X 2 + g 2 L6X 3 + g 2 5X 4 . (33) 

Supposing g 2 L 2 » 1, so that off-diagonal modes are heavy, we can ask when the two loop 
contribution to the partition function is parametrically smaller than the one loop partition 
function. If we normalise the one loop contribution to 1, then it is easy to see that the two 
loop contribution is of order g 2 N/(g 2 L 2 ) 2 . Thus the higher loop contribution is negligible 
if 

This condition is indeed satisfied by the width of the hemisphere distribution (9). Em- 
boldened, we go ahead and compute the one loop action for the simultaneous eigenvalues 
to obtain 

i i^j 

with 

V(s) = llog( s 2 )-ilog(l + 5 V). (36) 

In integrating out the off-diagonal modes one should gauge fix, for instance as described in 
[20]. This is of course the same action as we obtained before, except that now we have two 
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component vectors x = (x,y). The large N equations of motion are 

x= I l^V'{\x-y\)p{y)d 2 y (37) 
J y\ 

Writing this equation out for a radially symmetric distribution 

f L , . . f 2lT ,„ V'Ux 2 + r 2 - 2xr cos 0)(x - r cos 9) . . 

jo Jo V 2- ~\~ t — 2xr cos f 



X 

Jo Jo 

,2/„2 ^2 




! /r(r " = • ff(s r) I . (39) 



v / l + 2 5 2 (r 2 + x 2 ) + 5 4 i 



2 



In the second step we assumed, without loss of generality, that x > 0. We also introduced 

the step function H{s) which equals — 1 for s < and +1 otherwise. 

At large gL we can solve this last equation (39) with some educated guesswork. The 

solution is _ 

2N L 2 — r 2 

^) = 71 r- ( 4 °) 

where 

L = N'/ 2 (^y\ (41) 

The correctness of this solution can be checked by performing the integral (39) and then 
taking the strong coupling limit. This clearly does not agree with the hemisphere distri- 
bution (30) and the width L has a different scaling with A than (9). It is interesting to 
see that the eigenvalue problem (37) acquires a logarithmic non-analyticity in the 't Hooft 
coupling in the strong coupling limit. It is not easy to check the solution (41) numerically 
to high accuracy, because of the logarithmic dependence, but it easy numerically to see the 
distribution become paraboloid and the width scale like 1/A 1 / 4 rather than 1/A 1 / 6 . 

The implication of this mismatch is that the one loop effective action for the eigenvalues 
(35) is insufficient to capture the eigenvalue dynamics, despite the fact that the width (9) 
satisfies the naive bound (34). We suspect that this occurs because the scaling of L ~ 1/A 1 / 6 
is fairly close to the limiting scaling 1/A 1 / 4 . It may be possible to identify the required higher 
loop corrections to the effective action and re-obtain the correct emergent geometry. Indeed 
the appearance of a logarithm of the coupling suggests a resummation should be done, but 
we leave this for future work. 

To summarise this section: we have shown that the two matrix model (1) becomes 
commuting at strong coupling and that there is an emergent two dimensional hemisphere 
geometry (30). We saw that locality of physics in this geometry, i.e. that the off diagonal 
modes are heavy, was captured by the observables (27) and (28). We found that describing 
the eigenvalue dynamics will require going beyond the one loop effective action. 
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3 Solvable models with more than two matrices 



3.1 Bosonic model: No geometry 

There is a generalisation of the two matrix model that can also be treated analytically. It is 
instructive to see how this model does not lead to an emergent geometry at strong coupling. 
Consider the k + 1 matrix model with SO(k) symmetry only 

Z = J VXVY 1 ■ ■ ■ VY k e-^ x ' 2 -^ trY - + ^ ^ M W (42) 

This integral is again quadratic in all the Y m s, which may therefore be integrated out 
exactly. The crucial simplification is the absence of interactions between the Y matrices. 
Diagonalising X, we obtain 

Z = Jdx 1 ... dx N e~ ^ ^ M*i-^) 2 -f E w logli+^-a^] (43) 

The large N saddle point of this integral is an eigenvalue distribution satisfying 

f I 1 ~ k k 
J [x-y (x - y)[l + g 2 (x - y) 2 } 

We have rearranged the terms a little. As before f p{x)dx = N . The second term is a 
repulsive force of the same form as we found previously in (5). The first term, however, 
is now an attractive force (for k > 1) that is stronger at long distances. Eigenvalues 
separated by (x — y) 2 = \/{g 2 {k — 1)) experience no net force. One could presumably solve 
this integral equation fully using the techniques in [26]. Once again, we shall look for a 
pedestrian approach using the strong coupling expansion. 

A good starting point to attack (42) analytically is the limit of a large number k of 
matrices. The attractive force is becoming more important in this limit, so we might expect 
the width L of the distribution becomes small. In fact, based on the force balance equation, 
the natural scale to expect is L 2 ~ l/g 2 k. Let us assume this scaling and see where it leads. 

If L 2 ~ l/g 2 k and k is becoming large, then (x — y) 2 g 2 < 4L 2 g 2 is becoming small. 
Therefore we can expand the last term in (44). Rearranging, we obtain 

x + kg 2 £(-i) V n A* - y) 2n+1 P (y)dy = [ . (45) 

The first few terms in the sum read 

(1 + k\)x - ^(x 3 + 3a 2 x) + ^(x 5 + 10a 2 x 3 + 5a 4 *) + ■ ■ ■ = / • (46) 

JM iv z J x — y 



p(y)dy ■ 



(44) 
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We have used the fact that p(y) will be even and introduced the notation for the moments 
of the distribution 

a m = ^| y m P (y)dy. (47) 

Note that each term in this expansion is suppressed by a factor of L 2 g 2 ~ 1/k compared to 
the previous terms. Therefore we can solve order by order in the large k limit. The leading 
order solution is clearly the semicircle 

2N 



p(x) = —VL 2 ~^ 2 , (48) 



with 



This expression implies L 2 g 2 <C 1 for all A, in the large k limit. Therefore this is the 
correct leading order solution at large k for all couplings, not just strong coupling. It is not 
surprising that upon integrating out a large number of matrices one obtains the semicircle 
distribution, corresponding to effectively Gaussian degrees of freedom. We have checked 
that (49) agrees excellently with numerical results. 

To second order, the solution may be found using standard matrix model techniques. It 
is straightforward to verify that the solution to second order is 



N r- -/ k\ 2 L A k\ 2 I? 

where the width L satisfies (to this order in large k) 



r-z r/ k\ 2 Lr k\ 2 L 2 ,\ , 



4N 2 -2NL 2 (l + k\) + 3kL 4 X 2 = 0. (51) 

If we work in the large A limit for simplicity, this is seen to imply a corrected eigenvalue 
width of 

-/IK)- 

We checked these results by simulating the eigenvalue partition function (43) numerically. 
We again used Hybrid Monte-Carlo with N = 100 and stored 1000 configurations, deter- 
mining the distribution from iV x 1000 points. The resulting distribution for A = 1000 and 
k = 20 is shown in figure 3 together with the theoretical expectation (50). There is an 
excellent agreement. 

It is straightforward to move systematically to arbitrary order. For instance the seventh 
order solution is of the form 

p{x) = — ^Jh 2 - x 2 (a + bx 2 + cx 4 + dx 6 + ex 8 + fx 10 + gx 12 ) , (53) 

7T 
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Figure 3: Numerically simulated eigenvalue distribution solution to (43) with N = 100, 
k = 20 and A = 1000, together with the theoretical result to first (48) and second (50) order 
in a 1/k expansion. The leading order distribution has L = 0.0975 and the subleading 
distribution L = 0.1075. 

with the constants a, ■ ■ ■ ,g determined by plugging this expression into (45). The width of 
the distribution (at large A) is found to be at this order 

IW ( 3 1 5 1 111 57 1 33 1 39 1 \ 
~ V Afc V + 2fc + 4^ + TP + 32F + 64F + 16fc6y ' { ' 

For k = 2 this formula gives L = 2.3999-y/ N/X with small corrections. Given that k = 
2 is the smallest case we have to consider, this implies that (53) and (54) are a good 
approximation to the solution for all integer k > 1. The solution has L 2 ~ l/g 2 k, and 
so selfconsistently satisfies our scaling assumption. It seems plausible that the radius of 
convergence of the series (54) will be k = 1, consistent with the fact that we found a 
different scaling with A in that case. 

The implication of the above results, e.g. (54), is that for these multi-matrix models, 
with k > 2, the mass of the off diagonal modes of the matrix X is of the same order as the 
diagonal modes. This is because L 2 g 2 is order one or smaller. Therefore we do not expect 
an effective description in terms of the simultaneous eigenvalues of X and Y. Furthermore, 
we don't expect the Ys to commute amongst themselves as they only couple to one another 
through X. Let us see if these expectations are reflected in the observables we considered 
previously. 
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For simplicity, let us work to leading order at large A and k. That is, we take the solution 
(48) for p(x). We have seen that the scaling of quantities with A does not change away from 
this limit down to k = 2. 

The two point function of X is easy to evaluate 

tlX2 1 f L 2 / U L% N fx^ 

= — / x ptx)dx = — = — — . 55 

N N J_ L 4 2kX v ' 

To compute correlation functions involving the Ys, we introduce source terms as before 

Z[J}= J VXVY 1 ■ ■ ■ VY k e~ trX2 -^ trY ^ ^ m ^[x,Y m ] 2 +J2 m tr.j m Y m _ ^ 

Using exactly the same steps as in section 2.3 above, we can obtain the two point function 
for the Ys 

trY 2 1 f p(x)p(y)dxdy _ 1 f N 

-jr = Nj 2{l+ g^x-y) 2 ) * 2iV / P^P^ dxd y = 2 ■ < 57 > 

In this and subsequent formulae we are using the fact noted above that (x—y) 2 g 2 < 4L 2 g 2 <C 
1 for this distribution. We see that the Ys are spread out by an extra factor of A compared 
to the X matrix. The four point functions trX 2 Y 2 and tr XYXY are 



frX 2 Y 2 I p(x)p(y)x 2 dxdy 

- " J 2(1 • yM.r y) 2 ) 

L 2 N 2 N 3 



J 2(1 

\ j p( x )p(v) 



x 2 dxdy = = ^ , (58) 



and 



trXY n XY n 



I 



p(x)p(y)xydxdy 
2(l+g*{x-yy) 



This second expression is suppressed by an extra factor of 1/k and so can be neglected to 
leading order at large k. It follows that to leading order tr[X, Y] 2 = —2tvX 2 Y 2 . 
The four point functions involving two different Ys are 

-2 V 2 / p{x)p{y)p{z)dxdydz _1 f ^ ^ ^ Jm _ N 3 



*Y m Y n = J A{l+g 2 {x _ z?){l+g 2 {y _ zY) - 4 J P(x)p(y)p(z)dxdydz = — , (60) 
and 

N 

tiY m Y n Y m Y n = -~0. (61) 

This last term is sub leading in 1/N and therefore effectively zero in the planar approximation 
we are taking. The following four point functions are also obviously zero at large N, from 
the absence of Y mixing in the classical action 

tiXY m XY n = tiX 2 Y m Y n = . (62) 
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We have checked all the four point functions in this section by comparing with numerical 
results obtained by simulating the full partition function (42). 
We can now compute our ratios to be 

Ntx[X,Y n ] 2 = Ntv[Y m ,Y n } 2 

tiXHrY 2 trl^tri; 2 ' ^ 

and 

tr XY n XY n 1 
trX^ =k^ - (64) 

Both of these expressions are consistent with our expectation that this model is not com- 
muting. The fact that the first ratio is not going to zero and the second is less than one will 
remain true at finite k. It is interesting to note that at finite k the second ratio is not driven 
to zero as the coupling goes to infinity. Therefore the finite k models are not 'maximally 
non-commuting'. For the Ys however, we do have at large N that 

trYm Y n Y m Y n 



trYiY* 



0. (65) 



A simple but important lesson to draw from this model is that quantum loop effects are 
crucial in determining whether the strongly coupled system is commuting or not. Here, the 
extra one loop contributions from the matrices destabilised the emergent geometry of the 
two matrix model. This suggests that we can improve the situation by adding fermionic 
matrices to cancel the undesired loop contributions. 

A final point to note, from e.g. (46), is that the mass term is unimportant in deter- 
mining the strong coupling eigenvalue distribution in this model. This is in contrast to 
the (commuting at strong coupling) two matrix case we considered previously. Without a 
mass term, the coupling A can simply be absorbed into the normalisation of the X matrix 
and therefore does not have any dynamics associated to it. In particular, strong coupling 
cannot drive us to commutativity in a massless theory. A lesson to draw, therefore, is that 
commutativity should involve an interplay been the mass term and the 'commutator square' 
interaction terms. 

3.2 Model with fermionic matrices: still no geometry 

It has long been appreciated that supersymmetry facilitates the emergence of geometry [27, 
28, 29]. Here we show how a simple implementation of this idea works for us. However, the 
absence of direct interactions between the Y matrices will ultimately prevent the emergence 
of a geometry involving all the matrices in these solvable models. 
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We can supplement the 1 + k matrix model (42) with 2h fermionic fields as follows 
Z = J VXVY 1 ■ ■ ■ VY k VXi ■ ■ ■ V\ h Vm ■ ■ ■ Vfi h 

e -trX 2 -E m try^-E„trA nMn + 9 2 E m tr[X,y m ] 2 +i S E„trA n [X, Mn ] _ (gg) 

In this model, the X n and \x n are Hermitian N by N matrices of anticommuting numbers. 
Each matrix component of A and \x has only a single component, no 'spinor index' is 
necessary. Because of the asymmetry of the interactions, we are able to add fermions 
without introducing a Clifford algebra and without breaking any of the symmetries of the 
bosonic action. 

The partition function is quadratic in all the Ys and fermions, which we integrate out 
exactly to give 

Z = J dx! . . . dx N e- ^ M**-*i) 2 +*f* £*y logliWixi-xif] . (67) 

We see that including the fermions can cancel out the attractive part of the potential, thus 
increasing the likelihood that the eigenvalues of X will spread out sufficiently to give the off 
diagonal modes with the Ys a large mass. We already know the answer for the eigenvalue 
distribution in several cases: 

• If h < k — 1, the partition function is (42) with k' = k — h > 1 and one obtains the 
same results as in the previous subsection, with no emergent geometry. 

• If h = k — 1, the partition function becomes precisely that of section 2, leading to a 
parabolic distribution of the eigenvalues of X with width L ~ g^ 1 ^ 3 at strong coupling. 

• For the case h = k, the partition function is that for a Gaussian model and clearly 
results in a semi-circle distribution for the eigenvalues of X with width L ~ iV 1 / 2 . 

In the latter two cases, g 2 L 2 becomes large at strong coupling and therefore the off diagonal 
modes between X and any of the Ys become heavy. Thus they are candidates for emergent 
geometry. However, we will shortly see that the Ys will not commute amongst themselves. 
This is essentially an artifact of the simple action we have taken, with no interactions 
between the ys. Generically, if two Hermitian matrices Y m and Y n both commute with 
X then they should commute amongst themselves. The exception is if X has degenerate 
eigenvalues. This is effectively what is happening here. 

For the remaining cases, with h > k, all the terms in the effective action lead to repulsive 
forces. It is useful to write the large iV equations of motion as 



/ 



l+h-k ^ k-h 



x-y (x-y)[l + g 2 (x-y) 2 ] 



p(y)dy. (68) 
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Because the force is more repulsive than the case h = k, we expect that the distribution p(x) 
will have width L > N 1 / 2 . This would imply that gL S> 1 at strong coupling. Assuming 
this, then from (6) we can ignore the last term in (68). It follows that the distribution at 
strong coupling is a semicircle 

p( x )=™ v / U^, (69) 

7T.L 

with width 

L = y/2N(l + h-k). (70) 

This result includes the case h = k. The result self-consistently satisfies gL » 1 in the 
strong 't Hooft coupling limit. Therefore we can conclude that for all h > k — 1, the off 
diagonal modes of the Y matrices relative to X become heavy. 
It follows that 

trX 2 L 2 N(l + h- k) 

— = T = —2 (71) 

and, to leading order at large A using for instance (6) to compute the large A limit, 

tr*£ = _87V_ = 4V2N 1 
N 3vr 5 L 3vr(l + h- k) 1 / 2 X 1 / 2 ' 1 ' 

The commutator can be computed robustly as we did for the two matrix case, by differen- 
tiating the logarithm of the partition function with respect to A. To pick out the bosonic 
commutator squared, one should take the couplings in front of the bosonic and fermionic 
interactions to be distinct before differentiating. We find 

/V 3 

tr[X,y n ] 2 = -— . (73) 

Also proceeding as previously, we find the four point functions 

trX 2 Y 2 8LN 2 N 2 N 2 f3 
n ~ 15tt 5 Ag 2 + g 2 

8yg(l + h- kf^N 3 1 iV 3 N^f5_ 

15^ AV^"4A +_ T' ( } 



IhTig 2g z g l 
8y/2(l + h-k) 1 / 2 N 3 1 iV 3 

15^ AV^"2A +_ T' ( } 

2 2 3N 3 _ 3N 3 1 

trmn " W"8(lH-fc)A' (76) 
N 

txY m Y n Y m Y n = T ~0. (77) 
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As before, there is an unknown subleading contribution N 3 /3/\ due to the leading correction 
to the eigenvalue distribution at large A. We know that the contribution is equal in both 
the expressions above, because the difference has to reproduce the commutator (73). We 
do not need the value of (3 for our computations below, and we do not attempt to calculate 
it. 

The above results hold for h > k. In the 'critical' case h = k—1, the four point correlator 
of the Ys was not computed previously. Using the distribution (8) we find to leading order 
at strong coupling 

tr y 2 y 2 = 2 ^ 3 = ^ (*H J_ (m 
m B 280 5 2 L 2 140 \ 2 ) A 2 / 3 ' { ' 

We can now discuss the emergence of geometry in the models with h > k — 1, that 

is, with sufficiently many fermions that the off diagonal modes of X with any of the Y m 

matrices are heavy at strong coupling A — > oo. Two comments apply to all of these cases. 

Firstly, as we should expect, the X matrix commutes with the Y m matrices according to 

both of our criteria. Namely 

Nti[X,Y m } 2 txXY m XY m 
trXHrY* ~~ trX 2 ^ 2 ' aS (?9) 

This is easily seen from our above expressions for the relevant two and four point functions. 

Secondly, although one might therefore have expected the Ys to commute amongst 
themselves, this is not the case. For all of these models it follows from our above expressions 
that 

y ^W. ^y^ -O. as ( 8 °) 
11 1 m l n m n 

The fact that the Ys do not commute amongst themselves can be understood as being 
due to the absence of interactions between the Y matrices in the action of (66). At leading 
order in large N they are simply uncorrelated and cannot commute. Unfortunately, the very 
simplification that allowed us to get an analytic handle on this model undoes the possibility 
of an emergent geometry in which all of the matrices participate. The lesson we might take 
away is that genericity is another important property in the search for commutating models: 
if a matrix has degenerate eigenvalues it can commute with two other matrices without the 
other matrices needing to commute amongst themselves. 

Given that the X matrix commutes with all of the Y matrices, one might imagine 
picking one of the Y matrices, say Y\ (or perhaps some combination of them, to preserve 
the SO(k) invariance) and considering a geometry given by the joint eigenvalue distribution 
of this matrix with X. However, this geometry is not especially useful, even in the critical 
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case h = k — 1, where the X and Y\ eigenvalues would have roughly the same spread. 2 
The problem is that the physics of the remaining Y matrices would not be local in the y\ 
direction. Off diagonal modes of Y2 (say) connecting different values of y\ would have a 
comparable effect on the dynamics as the eigenvalues of Y2. So the strongly coupled matrix 
model is not solved by local two dimensional geometric physics in these cases. For emergent 
local geometry, all modes relating far away spacetime points should be massive compared 
to local modes. 

The natural next step, in search for higher dimensional geometry, is to introduce in- 
teractions between the Ys. This substantially increases the difficulty of solving the model 
analytically. 



4 The fully interacting multi-matrix model 

4.1 Bosonic model 

The interacting p-matrix bosonic model with full SO (p) invariance is 

Z = J VX 1 ■ ■ ■OT p e~ E " ,trI "' + ^ 2Em, ° tr[I "" X " 12 • ( 81 ) 

We cannot solve this model exactly, for p > 2, even in a strong coupling expansion. However, 
building on the intuition from previous cases, a few observations are possible. 

We do not expect this bosonic model to be commuting. We saw above that multiple 
one loop bosonic contributions result in an attractive potential at long distances. The 
eigenvalues were therefore insufficiently spread out for the off-diagonal modes to become 
parametrically heavier than the eigenvalues at strong coupling. Thus we do not expect 
an emergent geometry. We will study a commuting ansatz for a similar matrix model 
shortly. Let us suppose for the moment that all the elements of the matrices, diagonal and 
off diagonal, are of the same order. By SO(p) invariance, the entries will be of the same 
magnitude for all the matrices 

Here m c fj. is the effective mass of the matrix elements. If the matrices become heavy at 
strong 't Hooft coupling, then we might hope to trust a one loop evaluation of their masses 



m 2 eS . ~ I + g 2 y^iXikXk,} ~ -A- • (83) 

k m eff. 



2 The cases with h > k are even worse. The spread of the Yi eigenvalues is smaller by a power of A 
compared to the spread of the X eigenvalues, see (71) and (72). In commuting matrix models, a large 
anisotropy can prevent the smaller dimension from emerging at all [30, 31]. 
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From which we would conclude that the typical matrix entry scales as 

( 84 ) 

Note that if we diagonalise one of the X matrices, with entries of order (84), we will obtain 
diagonal elements of order A fl / 2 /A 1 / 4 . This is seen, for instance, from the observation that 
(84) implies trX 2 = XijXji ~ N 2 /A 1 / 2 whereas for a diagonal matrix trX 2 = X 2 ^ Thus 
(84) is precisely on the boundary of our condition (34) for trusting a one loop effective 
action. We will shortly give an independent argument supporting the scaling (84). 

In the previous paragraph, using selfconsistency and a non-commutativity assumption, 
we obtained a scaling for the spread of matrix elements with A at strong coupling. The 
SO(p) invariance of the model was also important; in section 3.1 above an anisotropic model 
gave different scalings. 

A consistency check of the above picture is that we can compute one of our ratios, 
following e.g. [15]. Consider the change of variables X m — > (1 + e)X m . This must leave 
the partition function invariant. To linear order the measure for each matrix changes as 
VX m — > (1 + (iV 2 — l)e)VX m . Requiring the partition function to be invariant and using 
the SO{p) symmetry leads to 

p(N 2 - 1) = 2p(trAl> - 2g 2 p(p - l)(tr[X m , X n } 2 ) . (85) 

This is an exact expression. Now using the scaling (84) at large coupling and taking the 
large iV limit we find 

trA^trX 2 2(p-l)g 2 (trX 2 J 2 ^ U ^ 
Thus, consistently, the matrices are indeed not commuting. 

At this point we can note a general result for models with these types of interactions. 
The model will be commuting, according to the observable (86), if and only if the eigenvalues 
are more spread out than N 1 / 2 \~ 1 / 4 as A — ► oo. Indeed, so far the only model we have 
discussed that satisfied this property for all the matrices involved was the two matrix model. 
The property is precisely the same condition that we found for the one loop effective action 
to be reliable. Although, as we saw in the two matrix case, this latter condition is not 
precise. 

Consistently with our discussion at the end of section 3.1 above for non-commuting 
models, we see that the mass term is not playing a role in any of our strong coupling 
considerations. See for instance (83) or (85). This allows us to make contact with previous 
numerical results on the model (81) in the absence of a mass term [12, 13, 15]. 
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If we naively take the strong coupling limit of the action (81), then we might expect 
to be able to drop the mass term. In the massless theory we could eliminate the coupling 
by rescaling the matrices X — > XA _1//4 . It is then immediate that the entries of X will 
scale like Xij ~ A _1//4 and therefore that eigenvalues will scale like x ~ 

A Tl/2 A -l/4_ Xhig 

scaling was indeed found numerically in [15] for p > 2. These results support our discussion 
above, including the observation that the mass term appears to indeed be unimportant at 
strong coupling in these models. For p = 2 the mass term is important because there is an 
infrared divergence in the massless model [12, 13] due to the zero modes in the commutator 
squared potential. This is consistent with the fact that we found a different scaling in the 
two matrix model case. 

The bottom line for these fully interacting bosonic models with p > 2 appears to be 
that the strong coupling physics essentially reduces us to the massless case. The scaling of 
quantities with A is fixed by dimensional analysis and the model is not commuting. Diagonal 
and off-diagonal modes contribute equally to generic observables, so there is no emergent 
local geometry. Clearly the arguments in this subsection are not intended to be as rigorous 
as our previous considerations. 

4.2 'Supersymmetrised' models 

We found in a solvable model above that adding fermions to the model such that the 
fermion determinant cancelled the bosonic one loop contribution allowed the eigenvalues to 
spread out. Supersymmetrisations of the massless fully interacting bosonic model exist in 
dimensions p = 3,4, 6, 10, corresponding to N = 2, 4, 8, 16 real supercharges, respectively. 
These have been studied numerically in for instance [13, 16]. In the cases of p = 6 and 
p = 10, the sign problem of the fermion determinants means that thinking in terms of a 
positive eigenvalue distribution is potentially misleading. However, in the cases of p = 3 
and p = 4 there is no sign problem at large N [11, 16]. Furthermore, the (massless) cases 
p = 3 and p = 4 show strong infrared effects: the partition function is divergent for p = 3 
whereas for p = 4 the partition function is finite but all moments (tvX 2m ) diverge [13]. 
Therefore these models are excellent candidates for an emergent commuting geometry upon 
adding a mass term, as we would expect a balance between the repulsive interactions and 
the confining mass term. 

Adding a supersymmetric mass term to the massless models is not straightforward, 
however. Instead we will consider simpler models. We start with the following SO(p) 
invariant matrix model, defined for all p, with 2(p — 1) Hermitian fermionic matrices A m 
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and fi m 

Z = j VX 1 ---VX p V\ 1 ---V\ p -iDn 1 ---Vn p - l 

- E m trX^ + I 3 2 E m ,„ tr[X m ,X„] 2 -E m A m ^l+ig 2 £JX„, .] 2 ^ m _ ^ 

The fermionic interaction has been chosen so that the bosonic and fermionic one loop 
determinants cancel. Note that 2(p — 1) is the appropriate number of massive fermions to 
cancel the massive bosons, because X m — ► UX m U~ l is a massless mode. In integrating out 
the off diagonal modes, one should for instance gauge fix as described in [20, 24]. 

Let us assume that this model is commuting at strong coupling and see if we can 
find a self consistent solution. If the off diagonal modes are heavy about this commuting 
background, we can hope that a one loop integration is sufficient. The eigenvalues are now 
vectors in W given by x\ = (xj, ■ ■ ■ , xf). The partition function becomes 

Z = J dx i . . . dx N e~ ^ 3 lo s l*-*il a . (88) 

The log \ xi — Xj\ 2 term can be thought of as a generalised Vandermonde determinant arising 
from simultaneously diagonalising the p matrices. The large N equations of motion are 



These equations of motion have been considered in some detail in, for instance, [20, 24, 31]. 
The solution is given by a (p — 1) sphere with constant eigenvalue density 



= bT-wffl- n a - t ' ■ (90) 



From the spherical eigenvalue distribution (90) one immediately obtains 

ixXl N 



(91) 



N 2p 

At strong coupling, the width of the eigenvalue distribution (90) is well within our bound 
(34) for when the one loop effective potential is expected to be valid. Therefore it seems 
likely that we have found a consistent saddle point for the full model. Without having 
solved the model exactly we cannot prove that this commuting saddle is the dominant large 
iV saddle of the integral. However, we can see that (89) describes a balance between a mass 
term and an eigenvalue repulsion of the type we found previously in a commuting model. 
The eigenvalues have spread out as far as the mass term allows and simultaneously made the 
off diagonal modes parametrically heavy at strong coupling. Therefore, there would seem 
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to be a good chance that this model is indeed showing an emergent geometry at strong 
coupling. 

The eigenvalue distribution (90) together with the validity of the one loop effective 
action allows us to compute the commutator 

1 f\y-y'\ 2 p(y)p(y')dPydPy' 



tr[X m ,X n ] 2 — f 

P J 



^ + g 2 \y-y'\ 2 

N 3 Xo\SP- 2 r (sinfl)P- 2 (l - cos#) 



df) 

pVo\SP~ l J l + A(l-cos0) 
N 3 1 

-> -, as A^oo. (92) 

p A 

Here we used the fact that the sphere has radius r 2 = N/2. We can also compute the 
leading order four point functions 

trX 2 m X 2 n = trX m X n X m X n = J # 'y p(y)y 2 m y 2 n = - = . (93) 

This computation simply involves integrations over the (p— 1) sphere. It should be possible 
to test these commuting saddle results using a numerical simulation of the full partition 
function (87). 

Other models with SO(p) invariance that have an a priori chance of having commuting 
saddles are the massless supersymmetric models we described at the start of this section 
together with a mass term for the bosons only: 

We can use the gamma matrix conventions of for instance [11]. In any case all that is 
important for us is that in the massless case the one loop bosonic and fermionic determinants 
cancel about a commuting background. We ignore for the moment the fermion zero modes. 
Therefore in our model, the one loop effective action about a commuting background is 

Z = jdx 1 . . .dx N e-^ +E ^ E^[iog|^-^l 2 -iog(i +3 2 |^-^| 2 )] . (95) 

The equation of motion for a radially symmetric eigenvalue distribution is therefore (p > 3) 

rL rn sm p - 2 6V'(\/x 2 + r 2 - 2xrcos6)(x - rcos9) 

da . 

V x 2 + r 2 — 2xr cos 9 

(96) 

In this expression the potential V is precisely as in (36) above. As in section 2.3 above, we 
can analytically perform the angular integral in (96) for any given p. 

We have not found an analytic solution for p(r) in (96). It is straightforward to see 
that if one inserts a generic ansatz for p(r) into the integral equation then the width of the 



[p - l)Vol^ 2 f drpir)^- 1 f 
Jo Jo 
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corresponding distribution will be L ~ 1/A 1 / 4 . This is on the borderline for selfconsistency 
of the one-loop action. Furthermore, we have simulated the commuting matrix model (95) 
to rather high accuracy in N and A. The results for the width of the eigenvalue distribution 
are shown in figure 4 for the cases of most interest, p = 3 and p = 4. 



N=500 i i 1 

N=1000 ^-x-^ 
N=2000 . 
N=3000 ' □ : 
-0.477log(>.)-0.277 . 



N=1000 h 
N=2000 - 
N=3000 r 
-0.488log(X)-0.175 



Figure 4: For the p = 3 model (left) an the p = 4 model (right), plots of log(A) against 
l°g(Sm J2iLi( x m) 2 /N 2 )- The fitting line is obtained from last two points, A = 8000 and 
A = 12000, with N = 3000. 



The results of figure 4 show that the width of the eigenvalue distribution at strong 
coupling is a little greater than L ~ 1/A 1 / 4 . This is reminiscent of what we found previously 
for the commuting bosonic two matrix model in section 2.3, where the width was enhanced 
by a log A factor. One possibility is therefore that these models are indeed commuting and 
that the eigenvalue dynamics is not captured by the one loop action (95). It would be 
interesting to test this possibility by simulating the full partition function (94) numerically 
and computing our observables (27) and (28). We cannot say much beyond this for this 
model with the results we have here. 

A curious observation we can make is that if we combine the two plots of figure 4 into 
one plot, then the two lines essentially lie on top of each other. 

One new feature that arises in this last model is that there are fermion zero modes 
about the commuting saddle. We did not add a mass term for the fermions and therefore 
fermionic matrices that simultaneously commute with the bosonic matrices have zero action. 
Therefore the partition function will vanish unless we insert some fermions to mop up the 
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zero modes. Specifically we should consider an insertion like 



<n*n---*iVAr>- 



(97) 



a 



We have not written this observable in an explicitly basis independent way, but rather in 
the basis given by the commuting Xs. The ^ matrices are anticommuting, so the observable 
has to be antisymmetrisable in the components of ^, as well as having N entries for each 
matrix. 

5 Classical fuzzy spheres 

There is another notion of geometry that appears repeatedly in studies of matrix theory and 
string theory. These are non-commutative spaces formed by the classical matrix degrees 
of freedom of collections of D branes. For instance, the DO brane matrices can form a 
regularized version of a higher dimensional membrane geometry (this was a crucial insight 
to develop matrix theory [27]). Such blowing up into higher dimensions often occurs due 
to the Myers effect [32], and it can also show up in the classification of vacua of some 
supersymmetric field theories [33]. 

A typical example of such spaces are fuzzy spheres, whose classical coordinates satisfy 
the SU (2) Lie algebra relations. There are many other such configurations that have been 
studied. In this section we will make a small study of fuzzy spheres, to complement our 
points of view from the rest of the paper. 

The first thing we would like is a matrix model whose classical saddle points give rise 
to fuzzy spheres. Define the following matrices, for i,j, k taking values 1,2,3: 



This action can be promoted to a potential for a subsector in a supersymmetric matrix 
quantum mechanical system (see for example the eleven dimensional plane wave matrix 
model [34]). A different approach to obtaining an emergent fuzzy sphere from a matrix 
model may be found, for instance, in [35, 36]. 

If we set m = 0, classically we have an infinite family of saddles forming a continuous 
space of solutions. These are characterized by = 0. This can be interpreted as a 



M k = ie ijk XiXj - mX k , 



(98) 



where the X are hermitian matrices. Now consider the action 



S = trM 2 



(99) 
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classical infrared divergence: the set of classical vacua is a non-compact manifold. There 
will be zero modes in this system and in principle the classical dynamics can run away along 
these directions if one adds some time dependence. 

Turning on the parameter m is adding a mass term, plus some other effect linear in m. 
These effects lift the flat directions and one ends up with a finite set of discrete classical 
vacua. These vacua are characterized by matrices of order m. The general solution is 
described by = 0. This can be solved by 

X\ = mLi, X 2 = W1L2, X3 = mL^ , (100) 

for some angular momentum representation {Li} of the Lie algebra of SU (2). Let us choose 
the N dimensional irreducible representation (the spin s = (N — l)/2 representation). 

We find in this case that the eigenvalues of X are of order N. We can evaluate our 
criterion for commutativity as a function of N. This is, let us compute the ratios 



-Ntr[X ll X 2 } 2 ti(X 1 X 2 X 1 X 2 ) 
n= trX*trX* ' r2= tr(^|) ' (101) 



^ = 1772 =77^-TT- (102) 



The first of these ratios is easy to compute. We get that X2] = mX^, and using the 
symmetry between X±, X 2 , X3, we conclude that 

N 3 
trZf ~ s(s + 1) 

When we take N to be large this expression vanishes and hence one can argue that the 
matrices approximately commute. 

The ratio r 2 is slightly harder to evaluate. We use several identities. Firstly: 

tr[Li, L 2 } 2 = -JV ^** 1 ) = 2tr(LiL 2 LiL 2 ) - 2tr(L?L^) . (103) 

We also find by using symmetries between the X matrices that 

N [s(s + l)] 2 = trL 2 . = 3tvLf + 6tr(L 2 Ll) , (104) 

and it is easy to show that 

trL? = —3(1 + s)(l + 2s)(-l + 3s + 3s 2 ) . (105) 
15 

With these results at hand, we compute that 

s 2 + s - 2 5 
r2= s2 + s + l/2 =1 - 2s(s + l) + r (1 ° 6) 
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We see that r2 approaches one as we take N large, again suggesting this can be considered 
as an approximately commuting matrix model. 

In contrast to our previous examples, here the spread of the eigenvalues is of order N 
rather than y/N. This large size can be considered as a non-perturbative effect, because 
there are many saddles of the matrix model. The claim is that in the large N limit these 
configurations of large representations go to a smooth geometry, while the vacuum where 
X = is non-geometric (as we have seen before). 

At this stage, one would also like to understand how to reconcile our picture of large 
off-diagonal masses causing commutativity with the fuzzy sphere case. Since the gauge 
group of the configurations is completely broken, and not U(1) N as in the case of strict 
eigenvalue saddles, something else must replace this notion. The answer can be found in 
the work of Bigatti and Susskind [37]. The idea is that in the presence of some background 
non-commutativity, the modes that have large momentum become extended, with an ex- 
tension proportional to the momentum and the noncommutativity. It is clear that here 
the momentum should be replaced by the angular momentum of fluctuations. We should 
therefore find that the fuzzy spherical harmonics correspond to stretched segments whose 
size grows with the angular momentum. The mass of these fluctuations also ends up being 
proportional to this angular momentum. This spectrum has been computed for various 
types of fuzzy sphere configurations (see for example [38, 39, 40]). It is clear that the modes 
with very high angular momentum should be very extended and massive, and for the most 
part the picture will not look too different from the massive off-diagonal modes that we have 
discussing so far in the other quantum models that we have solved. This seems to suggest 
that within string theory and matrix models, different notions of emergent geometry may 
be a lot closer than would appear at first sight. 

6 Summary and discussion 

The basic questions underlying this paper are as follows: suppose we are given a multi- 
matrix model that we cannot solve exactly. Two questions we might ask are the following. 
Do the degrees of freedom describing the model reduce in a strong coupling expansion from 
order N 2 to order N? Are these eigenvalue excitations governed by a simple (e.g. one loop) 
effective action? 

We have found that the answer to these questions depends on the details of the model. 
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In the cases we could solve exactly, we found that the observables 

Ntr[X,Y] 2 tiXYXY 

IrA^ and teXSyT' (1 ° 7) 

captured the property of whether the model became commuting at strong coupling. The 
first of these observables tends to zero for all pairs of matrices in a commuting model 
whereas the second tends to one. We found that this occurred in a bosonic two matrix 
model with a commutator squared interaction. We showed that this model is described 
at strong coupling by an emergent two dimensional hemisphere geometry. The eigenvalue 
dynamics about this geometry, however, is not describable by the one loop effective action. 

For the other solvable models we considered there was no emergent geometry at strong 
coupling. In the bosonic model of section 3.1 this occurred because integrating out the off 
diagonal modes induced a strong attractive force on the eigenvalues. The clumped together 
eigenvalues then had light off-diagonal modes connecting them. In the fermionic model of 
section 3.2 this attractive force was cancelled, but the lack of interactions between some 
of the matrices (necessary to solve the model exactly) meant that these matrices were 
uncorrelated and did not commute. 

We then turned to fully interacting models which we cannot solve exactly. For the 
bosonic model we argued that there was no local emergent geometry. For models with 
fermions we found apparently consistent commuting saddles in a model in which the fermionic 
measure was constructed by hand to cancel the bosonic one loop effective action. These 
models appear to have an emergent spherical geometry. It would be very interesting to test 
our (self-consistent) results for that model with a full numerical simulation of the partition 
function (87). For the supersymmetric model with a bosonic mass term (94) we did not 
find clearcut evidence for emergent geometry. It would also be very interesting to simulate 
that model numerically and to compute our observables (107) at strong coupling. 

We also considered some classical fuzzy spaces as other toy models where geometry can 
be argued to appear. We found that in these cases, at the classical level one also found that 
the matrices making up the fuzzy sphere are close to being commuting in the large N limit. 
The size of the matrices ends up being of order N rather than y/N \ and this is what makes 
them more classical. We also found that in these systems the notion of off-diagonal modes 
should be captured by the spherical harmonics of the fuzzy sphere. The angular momentum 
should be correlated with the size of the segment joining two putative points on the sphere. 

An important feature of all the multi-matrix models (both quantum and classical) that 
we have considered in this paper is that they have a mass term. This is very natural in 
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AdS/CFT, because the boundary of AdS in global coordinates contains a spatial three- 
sphere which gives even the zero modes a conformal mass. It is different, however, from 
previous considerations of matrix models as putative nonperturbative formulations of string 
theory, as these do not usually have mass terms. In the examples we considered, a balance 
between the mass term and the commutator squared interaction term was important for 
the emergence of local commuting geometry at strong coupling. Otherwise the coupling 
dependence is simply dimensional analysis. In particular, the cases with emergent geometry 
occur when the massless model has infrared divergences. These infrared divergences signal 
a potential instability of the system to grow in size. It is the regulation of this growth that 
seems to give us a notion of geometry. If the size is not controlled, one can imagine that 
the end of such a scenario is a system where all the eigenvalues have scattered at infinite 
distance from each other and there is nothing interesting left. We would definitely not want 
to call such type of configuration a geometry, but the dynamical process of reaching such 
an end configuration could be a cosmology of sorts. Previous cosmological applications of 
matrix models include [41, 42, 43]. 

We believe it is worth exploring these ideas further. There are various supersymmetric 
matrix quantum mechanical models where there is a mass term (for example, the plane wave 
matrix model [34]), and our notions of geometry might be useful to tackle the spectrum of 
extended objects in such systems. In a different direction, recent progress has been made 
connecting black hole physics with supersymmetric matrix models at finite temperature 
(with infrared divergences) [44, 45, 46]. Infrared divergences can also be found in weakly 
coupled N = 4 Super Yang-Mills theory on a three sphere, e.g. [47]. At high temperatures 
one expects off diagonal modes to play an important part in the dynamics. Therefore one 
faces an interesting challenge in combining our notion of an emergent geometry with black 
hole physics. 
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A Some corrections away from large A 
A.l Correction to the parabolic distribution 

In order to compute the four point functions in section 2.2 to subleading order, we needed 
the leading order correction to the large A eigenvalue distribution 

P(x) = Po(x) + ^ . (108) 

Here po(x) is the parabolic distribution (8). It turns out that for the leading order correction 
we can take L to remain given by (9) and f p\(x)dx = 0, thus retaining J p(x)dx = N. 

Using (6) and expanding the equations of motion (5) to second order in strong coupling, 
we find the following expression for the correction to the distribution 

Clearly this solution is somewhat formal, as the eigenvalue density diverges as x — ► ±L. 
However, inside the integral it does solve the expanded equation of motion. Furthermore, 
the region in which the full distribution (108) becomes large is exponential small (~ e _Al/3 ). 
This is related to the fact that the parabolic solution po(x) was not valid very close to the 
endpoints, as we noted in the main text. Therefore, we should be able to use (109) to 
evaluate observables for which the integral over the eigenvalues is finite. 

Using the integrals in (22) and (23) it is easy to use (109) to find the correction to the 
four point functions that is due to pi (x) : 

N 3 

5(trXYXY) = 5{trX 2 Y 2 ) = -— . (110) 
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